##########################################################################################

library(data.table)
library(optparse)
library(dplyr)
library(ggplot2)
library(parallel)

##########################################################################################

option_list <- list(
    make_option(c("--sample_list_file"), type = "character"),
    make_option(c("--diff_file"), type = "character"),
    make_option(c("--gtf_file"), type = "character"),
    make_option(c("--pathway_path"), type = "character"),
    make_option(c("--q_t"), type = "character"),
    make_option(c("--foldchange_t"), type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    sample_list_file <- "~/20220915_gastric_multiple/rna_combine/analysis/config/tumor_normal.list"
    diff_file <- "~/20220915_gastric_multiple/rna_combine/analysis/images/DiffGene/DiffGene.tsv"
    out_path <- "~/20220915_gastric_multiple/rna_combine/analysis/images/DiffGene"
    pathway_path <- "~/ref/Pathway/"
    gtf_file <- "~/ref/GTF/gencode.v19.ensg_genename.txt"
    q_t <- 0.05
    foldchange_t <- 2

}

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

sample_list_file <- opt$sample_list_file
diff_file <- opt$diff_file
out_path <- opt$out_path
pathway_path <- opt$pathway_path
gtf_file <- opt$gtf_file
q_t <- as.numeric(opt$q_t)
foldchange_t <- as.numeric(opt$foldchange_t)

dir.create(out_path , recursive = T)

##########################################################################################

info <- data.frame(fread(sample_list_file))
result_wilcox <- data.frame(fread(diff_file))
dat_gtf <- data.frame(fread(gtf_file , header = F))
colnames(dat_gtf) <- c("gene_id" , "Hugo_Symbol")
result_wilcox <- merge( result_wilcox , dat_gtf , by = "gene_id" )

##########################################################################################

class_type <- c( "Normal" , "IM" , "IGC" , "DGC")

##########################################################################################

fisher_result <- function(P0.05_all,bg_gene,gene_pathway,pathway_name){
    ta <- length(which(P0.05_all%in%gene_pathway))
    tc <- length(which(bg_gene%in%gene_pathway))
    tb <- length(P0.05_all) - ta
    td <- length(bg_gene) - tc
    data_fisher <- matrix(c(ta,tb,tc,td),nrow=2)
    fisher_res <- fisher.test(data_fisher)
    
    gene_P_in_pathway <- P0.05_all[which(P0.05_all%in%gene_pathway)]

    if(length(gene_P_in_pathway)==0){
        gene_paste <- ""
    }else if(length(gene_P_in_pathway)==1){
        gene_paste <- gene_P_in_pathway
    }else{
        gene_paste <- gene_P_in_pathway[1]
        for(i in 2:length(gene_P_in_pathway)){
            gene_paste <- paste0(gene_paste,",",gene_P_in_pathway[i])
        }
    }

    res <- data.frame(pathway=pathway_name,enrich_top=paste(ta,length(P0.05_all),sep="|"),enrich_all=paste(tc,length(bg_gene),sep="|"),bggene_in_pathway=tc,pathway_gene=length(gene_pathway),P=fisher_res$p,OR=fisher_res$estimate,gene=gene_paste)
    res
}


computPathway <- function(pathway_path = pathway_path , P0.05_all = P0.05_all , bg_gene = bg_gene , direction = direction , output = output ){
    
    dir.create(output , recursive = T)

    for( pathway_file in grep( "symbols.txt" , list.files(pathway_path) , value = T) ){

        file <- paste0(pathway_path , "/" , pathway_file)
        dat_pathway <- fread(file)

        pathway <- unique(dat_pathway$pathway)
        pathway_res <- data.frame(rbindlist(mclapply(pathway,function(x){
            # print(x)
            gene_pathway <- subset(dat_pathway,pathway==x)$gene
            #gene_pathway_procod <- gene_pathway[which(gene_pathway%in%bg_gene)]
        
            res_enrichment <- fisher_result(P0.05_all,bg_gene,gene_pathway,x)
            return(res_enrichment)
        },mc.cores=20)))
        results <- subset(pathway_res,OR>1)
        resultso <- results[order(results$P),]
        resultso$P_fdr <- p.adjust(resultso$P,method="fdr",n=nrow(resultso))

        out_file <- gsub("[.]symbols[.]txt" , "" , pathway_file)
        write.csv(resultso,paste0(output,"/",out_file,".",direction,".csv"),row.names=FALSE)
    }
}

##########################################################################################
## 分高表达和低表达，分布做通路富集
for( pathology in unique(result_wilcox$class1) ){
    print(pathology)

    use_dat <- subset( result_wilcox , class2=="Normal" & class1==pathology)

    highExpression_gene <- subset(use_dat , padj < q_t & log2FoldChange > log2(foldchange_t) )$Hugo_Symbol
    lowExpression_gene <- subset(use_dat , padj < q_t & log2FoldChange < log2(1/foldchange_t) )$Hugo_Symbol

    ## IM相比于Normal显著差异的基因
    type <- paste0(pathology,"_Normal")
    bg_gene <- use_dat$Hugo_Symbol
    output <- paste0(out_path , "/" , type)

    ## 上调基因
    P0.05_all <- highExpression_gene
    direction <- "up"
    computPathway(pathway_path = pathway_path , P0.05_all = P0.05_all , bg_gene = bg_gene , direction = direction , output = output )

    ## 下调基因
    P0.05_all <- lowExpression_gene
    direction <- "low"
    computPathway(pathway_path = pathway_path , P0.05_all = P0.05_all , bg_gene = bg_gene , direction = direction , output = output )
}

##########################################################################################
## 关注HallMarks通路
dat_plot <- c()
for( pathology in unique(result_wilcox$class1) ){
    print(pathology)

    use_dat <- subset( result_wilcox , class1=="Normal" & class2==pathology)
    type <- paste0(pathology,"_Normal")
    output <- paste0(out_path , "/" , type)

    up_hall_file <- paste0( output , "/h.all.v7.5.1.up.csv")
    down_hall_file <- paste0( output , "/h.all.v7.5.1.low.csv")

    dat1 <- data.frame(fread(up_hall_file))
    dat2 <- data.frame(fread(down_hall_file))
    if(nrow(dat1) > 0 & nrow(dat2) > 0){
        dat1$direction <- "Up"
        dat2$direction <- "Down"
        dat_tmp <- rbind(dat1 , dat2)
        dat_tmp$pathology <- pathology
        dat_plot <- rbind(dat_plot , dat_tmp)
    }   
    
}

## 显著通路
dat_plot <- subset( dat_plot , P_fdr < 0.05 )
dat_plot$pathway <- gsub( "HALLMARK_" , "" , dat_plot$pathway )
dat_plot$pathology <- factor( dat_plot$pathology , levels = c("IM" , "IGC" , "DGC") )

## 存在同一个通路在显著上升和下机的基因中均显著
dat_plot_2 <- c()
for( path in unique(dat_plot$pathway)){
    tmp <- subset( dat_plot , pathway==path )
    if(length( unique(tmp$direction) ) > 1){
        print(path)
        tmp$pathway <- paste0( tmp$pathway , "(" , tmp$direction , ")" )
    }

    dat_plot_2 <- rbind( dat_plot_2 , tmp )
}

up_pathway <- dat_plot_2[dat_plot_2$direction=="Up",]
up_pathway_order <- unique(up_pathway[order(up_pathway$P_fdr , decreasing = T),"pathway"])
low_pathway <- dat_plot_2[dat_plot_2$direction=="Down",]
low_pathway_order <- unique(low_pathway[order(low_pathway$P_fdr , decreasing = T),"pathway"])
pathway_order <- unique(c(low_pathway_order , up_pathway_order))
dat_plot_2$pathway <- factor( dat_plot_2$pathway , levels = pathway_order , order = T)

image_name <- paste0( out_path , "/DiffGene.Hallmarks.pdf" )

col_direction <- c("blue" , "red")
names(col_direction) <- c("Down" , "Up")

p <- ggplot(dat_plot_2, aes(x=pathology, y=pathway)) + 
  geom_point(aes(size=-1*log10(P_fdr),color=direction))+
  scale_color_manual(values = col_direction)+ 
  theme_bw()+  #设置背景
  labs(size="-log10(Q.value)") +
  theme(panel.background = element_blank(),#设置背影为白色#清除网格线
    legend.position ='right',
    panel.grid.major=element_line(colour=NA),
    plot.title = element_text(size = 8,color="black",face='bold'),
    legend.text = element_text(size = 10,color="black",face='bold'),
    axis.text.y = element_text(size = 10,color="black",face='bold'),
    axis.title.x = element_text(size = 10,color="black",face='bold'),
    axis.title.y = element_text(size = 10,color="black",face='bold'),
    axis.ticks.x = element_blank(),
    axis.text.x = element_text(size = 10,color="black",face='bold') ,
    axis.line = element_line(size = 0.5)) 
ggsave( image_name , p , height=10, width=6) 